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INVESTIGATION  OF  THE  APPLICATION  OF  "ARRAY 
ALGEBRA"  TO  TERRAIN  MODELING 

by 

James  R.  Jancaitls  and  Ronald  L.  Magee 

Mapping  Developments  Division 
Topographic  Developments  Laboratory 
US  Army  Engineer  Topographic  Laboratories 
Fort  Belvoir,  Virginia  22060 


Background. 

One  of  the  missions  of  the  Mapping  Developments  Division  of  the  US 

Army  Engineer  Topographic  Laboratories  (ETL)  is  to  investigate  new 

techniques  and  equipment  which  may  improve  the  Defense  Mapping  Agency's 

cartographic  capabilities.  This  Includes  the  development  of  new 

cartographic  products,  the  Improvement  of  product  quality,  increasing 

production  rates,  and  the  addition  of  new  capabilities  to  existing 

products.  One  such  product  now  under  development  at  ETL  is  an  improved 

automated  digital  terrain  modeling  procedure  based  upon  the  Weighting 

Function  Interpolation  Technique  (WIT)^,  and  which  involves  the  generation 

of  numerous  polynomials  calculated  subject  to  the  least-squares  criteria. 

WIT  is  an  Integral  part  of  the  software  entitled  Contouring  Via  the 
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Surface  Averaging  Concept,  or  CONSAC  , which  is  ETL's  digital  terrain 
modeling  and  automated  contouring  computer  program.  The  terrain  model 
constructed  by  CONSAC  Includes  a rigorous  mathematical  smoothing  of 
observation  errors  and  results  in  a compact  and  continuous  functional 
model  over  arbitrarily  large  areas.  The  contours  are  produced  by 


3 


analytical  solution  of  the  Intersection  of  constant  horizontal  planes 
with  the  functional  terrain  model  which  guarantees  smooth,  continuous, 
arid  appropriately  semi-parallel  lines.  The  key  to  the  efficiency 
and  applicability  of  this  software  is  the  sequential  application  of 
locally  valid,  simple  functional  forms.  The  sequential  nature  of  these 
algorithms  results  in  the  capability  to  consistently  model  arbitrarily 
large  regions.  The  use  of  simple  functional  forms  for  the  basic  terrain 
model  allows  for  the  efficient  derivation  of  desired  terrain  Information, 
such  as  the  contour  lines.  The  DMA  production  applications  impacted  by 
this  software  include: 

a.  Orthophotographs 

b.  Automated  contouring 

c.  Weapons  systems 

d.  Databases 

This  software  accomplishes  the  basic  smoothing  and  modeling  tasks 
by  the  utilization  of  least  squares  linear  approximation  techniques  on 

3 

the  gridded  elevation  data  produced  by  automatic  compilation  equipment. 

Introduction. 

4 

In  1972  Dr.  Rauhala  introduced  an  Improved  algorithm  for  linear 
approximation  of  data  lying  on  an  orthonormal  grid  subject  to  the  least- 
squares  criteria  and  subsequently  Included  it  in  his  1976  "Array  Algebra" 
manuscript.^  This  work  and  the  further  investigations  done  by  Kratky 
(1976)  and  Crombie  (1976  a and  b) ^ ’ came  to  our  attention  in  June  1976. 
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Kratky  stated  that  his  efforts  involving  polynomial  transformation  of 
satellite  imagery  enjoyed  a significant  increase  in  computational 
efficiency  and  accuracy  as  a result  of  the  inclusion  of  this  new 
modification  to  the  least-squares  algorithm.  Based  upon  the  pioneering 
work  of  Rauhala  and  the  subsequent  success  of  Kratky  and  Crombie,  this 
technique  was  considered  to  possibly  have  a cost-effective  application 
to  the  approximation  algorithm  in  ETL's  terrain  modeling  procedure. 

A direct  implementation  was  not  undertaken  Immediately  due  to  some 
uncertainties  we  had  in  the  theoretical  basis  and  applicability  of  this 

technique  to  our  problem.  At  that  time,  Rauhala' s proof  of  the  array 
algebra  algorithm's  equivalence  to  the  conventional  least-squares 
algorithm  was  difficult  to  follow  due  in  part  to  poor  English  language 
translations  of  Rauhala' s work,  and  in  part  to  his  deviation  from  the 
simpler  currently  accepted  matrix  notat^.onal  forms  in  those  cases  where 
they  could  have  been  utilized.  At  ETL  the  proof  of  these  algorithms^ 
equivalence  was  first  investigated  by  Crombie  (1976  a),  who  produced  a 
lengthy  derivation  of  their  algebraic  equivalence  using  a laborious 
presentation  and  manipulation  of  the  basic  system  of  scalar  equations. 

Of  concern  to  us  was  the  fact  that  the  derivation  of  the  array  algebra 
approximation  algorithm  did  not  proceed  from  an  a priori  stated 
mathematical  condition,  as  does  the  least-squares  normal  equations, 
splines,  and  all  other  well-used  approximation  algorithms.  Further, 
the  computational  efficiency  described  by  Kratky,  Rauhala,  and  Crombie 
(1976  b)  Involved  cases  similar  to,  but  not  exactly  like,  our  application. 
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For  the  above  reasons,  it  was  determined  that  the  application  of  array 
algebra  to  ETL's  terrain  modeling  procedure  was  to  be  Investigated  in 
the  following  manner: 

1.  A theoretical  analysis  of  array  algebra  would  be  performed  in 
order  to  verify  specifically  the  equivalence  of  array  algebra  and  the 
conventional  least-squares  solutions  for  our  application. 

2.  A paper  study  would  be  conducted  to  compare  the  computational 
efficiency  of  ETL's  terrain  modeling  algorithm  using  the  current 
conventional  least-squares  method  and  the  array  algebra  technique. 

3.  An  empirical  study  involving  the  development  of  two  digital 
computer  programs  would  be  conducted  to  verify  the  paper  study  results, 
and  to  evaluate  those  aspects  of  computer  implementation  not  easily 
predicted,  such  as  FORTRAN  "DO-LOOP"  overhead. 

4.  Given  that  array  algebra  proved  to  be  a cost-effective  technique 
for  ETL's  terrain  modeling  procedure,  it  would  then  be  implemented 

into  CON SAC. 

In  the  following  sections,  this  paper  reports  the  status  and  results 
of  these  investigations.  In  order  to  facilitate  the  presentation  of 
this  paper,  all  of  the  mathematical  derivations  have  been  placed  in  the 
appendixes  with  only  the  major  results  contained  in  the  bodv. 
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Theoretical  Analysis 

As  stated  in  the  introduction,  one  of  our  major  initial  concerns 
regarding  the  "array  algebra"  algorithm  was  its  mode  of  derivation. 


6 


In  contrast  to  all  other  widely  used  numerical  algorithms  , the  array  ' 

algebra  solution  did  not  proceed  from  an  a priori  stated  mathem.r  ical 

condition.  Rauhala's  derivation  of  the  solution  equations  for  the  | 

simplest,  two-dimensional  case,  involves  a utilization  of  the  bilinear  J 

form  in  concert  with  only  matrix  manipulations,  (see  appendix  C).  This  ^ 

method  of  derivation  is  unillumlnatlng  in  terms  of  what  this  solution 

1 

does  for  the  problem  solver  employing  it.  The  previously  mentioned 

proofs  of  equivalence  of  the  array  algebra  and  least-squares  solutions  i 

are  equally  unilluminating  derivations.  For  these  reasons,  this  effort 

was  directed  toward  a proof  and  verification  of  these  solutions 

equivalence  which  followed  more  traditional  lines.  More  specifically,  a 

derivation  of  the  two-dimensional,  array  algebra  solution  was  sought,  which 

was  a systematic  derivation  of  the  numerical  solution  that  is  easily 

seen  and  proceeds  from  an  a priori  stated  mathematical  condition  using 

direct  and  simple  analyses  based  upon  widely  known  and  understood  basic 

principles. 

Such  a derivation  was  found  anc  will  now  be  presented  in  the  following 
format.  On  the  left-hand  side  of  the  page  the  well-known  derivation  of 
the  least-square's  normal  equations^ ^ will  be  presented,  and  on  the 
right-hand  side, the  corresponding  steps  of  our  new  derivation  of  the  two- 
dimensional  array  algebra  solution  of  interest  in  our  application.  (For 
a complete  definition  of  all  matrices  see  appendixes  A and  C) . First, 
a model  equation  is  chosen, with  the  array  algebra's  bilinear  form 
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dependent  upon  the  observations  lying  on  an  orthonorraal  grid  in  some 
coordinate  system, 

Z = AC  W = XC 

L L A 

and  the  objective  function  formed 

4)  = (Zl  - AC^)T(Z  - AC^)  4)  = llw  -XC^Y^ll 

then  minimized  by  differentiation 

li  = (Z,  - AC,  )'^(-2A)  = X^(W  - XC.Y^)Y(-2) 

SCl  L L 3Ca  ^ 

setting  these  partials  equal  to  zero  and  solving  by  rearranging  and 
assuming  the  existence  of  inverses  gives 

T -1  T T -1  T T -1 

C,  = (A  A)  C.  = (X  X)  X^WY(Y  Y) 

u L A 

The  one-to-one  correspondence  of  steps  and  similarity  in  dependence 
upon  the  application  of  calculus  of  variations  results  in  a much  more 
illuminating  theoretical  basis  for  the  array  algebra  solution.  This  new 
derivation  also  represents  a much  clearer  proof  of  the  equivalence  of 
these  algorithms  because  a simple  expansion  of  the  model  and  objective 
functions  to  their  basic  scalar  equation  form  demonstrates  their  equality. 
As  is  well-known,  two  solutions  of  an  approximation  problem  are  equivalent 
provided  the  model  equations  and  objective  functions  are  equal.  The 
derivation  presented  here  also  allows  for  an  easy  and  straightforward 
Inclusion  of  "structured"  observation  weighting  and  arbitrary  control 
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constraints.  The  modification  of  the  algorithm  for  Kalman  Filter 
processing  has  not  yet  been  made  but  seems  to  be  a definite  and  useful 
possible  extension  at  this  time.  Appendix  A contains  a complete  and 
detailed  presentation  of  our  derivation  of  the  two-dimensional  array 
algebra  solution,  as  well  as  the  derivation  when  "structured"  observational 
weighting  and  observation  constraints  are  included.  Appendix  B contains 
the  detailed  derivation  of  the  standard  least-squares  normal  equations,  as 
well  as  a numerical  example.  Appendix  C contains  a detailed  presentation 
of  Rauhala's  derivation  of  the  two-dimensional  array  algebra  solution, 
as  well  as  the  same  numerical  example  of  appendix  B,  only  this  time  solved 
by  the  array  algebra  algorithm. 

In  conclusion,  the  two-dimensional  array  algebra  solution  of  Rauhala 
represents  another  computational  procedure  for  approximation  subject  to 
the  least-squares  critera  when  the  observations  lie  on  an  orthonormal 
grid  in  some  coordinate  system.  We  have  presented  a derivation  of  this 
algorithm  which  follows  more  traditional  lines  and  is  more  illuminating 
with  regard  to  what,  exactly,  the  array  algebra  algorithm  does  for  the 
problem  solver  utilizing  it, 

1 

Computational  Requirements  Analys is 

The  computational  requirements  analysis  was  completed  in  October  1976, 
and  provided  some  insight  as  to  the  Impact  array  algebra  would  have  on 
ETL's  terrain  modeling  algorithm,  CONSAC  II.  In  order  to  isolate  the 
computational  differences  between  the  conventional  least-squares  and 
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array  algebra  techniques  we  essentially  performed  a quantitative 
comparison  of  the  number  of  computer  multiplications  and  additions 
needed  to  compute  the  coefficient  matrices  or  where: 

T “IT 

= (A  A)  A Zj^,  the  least -squares  normal  equation 
or 

= (X^X)“^X^  Y(y'^Y)“^  , the  array  algebra  solution 

and:  Z,  and  Z,  are  matrices  of  observation 

L A 

A is  a matrix  of  polynomial  terms 
X and  Y are  uni-directional  parameter  matrices 

(for  further  specification  see  appendixes  B and  C) . 

Since  each  C,  or  C.  describes  only  one  polynomial,  then  it  would 
L A 

seem  that  it  must  be  calculated  once  for  every  polynomial  needed  to 

describe  an  entire  map  sheet.  Fortunately,  this  is  not  true.  Since 

(A  A)”  A^  is  dependent  only  upon  the  x-y  coordinates  of  each  of  the 

T -IT 

observations  being  fit,  then  (A  A)  A must  be  calculated  whenever  the 
x-y  location  of  the  set  of  observations  change.  If,  however,  we  choose 
to  move  the  origin  of  the  entire  coordinate  system  commensurately  with 
our  progression  through  the  set  of  observations,  then  (A  A)  A is 

constant  and  independent  of  the  position  of  the  origin.  Thus,  use  of 

T —IT 

this  local  roving  coordinate  system  allows  us  to  calculate  (A  A)~^A 
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only  once  per  map  sheet.  The  same  technique  prevails  for  (X^X)”^X^ 
and  Y(y'^Y)~^  . 

Now  it  can  be  seen  that  both  array  algebra  and  the  conventional 
least-squares  technique  require  two  computational  events:  one-time 
operations  and  repetitive  operations . (see  figure  1.) 


LEAST-SQUARES 

ARRAY  ALGEBRA 

ONE-TIME 

OPERATIONS 

M = (A^A)~^A^ 

N = (x'^X)"^\'^ 

T - 1 

Q = Y(Y  YJ 

REPETITIVE 

OPERATIONS 

C.  = NZ.Q 

A A ' 

Figure  1. 

One-Time  and  Repetitive 
Operations 

The  one-time  operation  is  the  calculation  of  the  solution  matricies 
in  the  local  roving  coordinate  system.  The  repetitive  operation  involves 
the  use  of  the  solution  matrices  repetitively  for  each  subarea  of 
the  terrain  being  modeled  locally. 

The  rigorous  derivation  of  the  multiplication  and  addition  times 
involved  in  these  operations  is  treated  in  appendix  D. 

A calculation  of  the  difference  of  the  number  of  one-time  multiplications 
required  by  the  conventional  least-squares  method  and  the  array-algebra 
method  produced  the  results  shown  in  figure  2a.  Notice  that  the 
computational  savings  Increase  with  the  size  of  the  model  polynomial  and 
with  the  number  of  elevation  observations  per  polynomial  fit.  One-time 
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V = The  difference  of  least-squares  one-time  multiplications  and  array 
algebra  one-time  multiplications. 

N = Number  of  observations  per  polynomial 

M = Number  of  terms  in  model  polynomial 


Figure  2a:  Difference  of  one-time  multiplications 


2 


addition  Is  similar  and  is  shovm  as  figure  2b. 

The  same  difference  concerning  total  repetitive  multiplications 
yielded  similar  results,  as  shovm  in  figure  3a.  Again,  repetitive 
additions  were  similar  and  are  shown  in  figure  3b. 

Given  that  the  modeling  procedure  is  run  on  a CDC-6400  computer 
and  limiting  the  number  of  elevation  points  per  polynomial  to  an  empirically 
determined  optimal  value  (see  appendix  E) , then  the  computer  multiplication 
and  addition  time  was  seen  to  vary  with  the  size  of  the  model  equation. 
Figure  4 demonstrates  the  computational  superiority  of  the  array  algebra 
technique  as  it  applies  to  the  CONSAC  II  terrain  modeling. 

This  study  shows  that  the  smaller  the  local  polynomials  the  more 
efficient  the  terrain  modeling  procedure.  Note  that  the  fact  that  more 
complex  polynomials  can  cover  larger  areas  to  within  a given  accuracy 
has  been  included  in  this  analysis.  Previous  research  had  demonstrated 
this  fact  and  the  current  ETL  terrain  modeling  software  makes  use  of  the 
most  efficient  four  coefficient  model  equation.  At  this  point  on  the 
curve,  the  savings  ratio  predicted  by  this  paper  study  is  1.75  to  1 if  the 
array  algebra  algorithm  is  employed  Instead  of  the  least-squares  normal 
equations. 

Empirical  Analysis. 

In  order  to  more  thoroughly  investigate  these  algorithms'  relative 

computation  requirements  an  empirical  study  was  undertaken.  This 

( 

empirical  study  was  deemed  necessary  to  fully  and  unambiguously  compare 
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V = The  difference  of  least-squares  one-time  additions  and  array 
algebra  one-time  additions. 

N * Number  of  observations  per  polynomial 

M = Number  of  terms  in  model  polynomial 

Figure  2b:  Difference  of  one-time  additions 
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Figure  3b:  Difference  of  Least-Squares  and 
Array  Algebra  repetitive  additions  per  map  sheet 
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Figure  A:  Theoretical  comparison  of  efficiency: 
Least-Squares  vs  Array  Algebra 
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those  aspects  of  these  algorithms  which  are  not  as  easily  studied 
analytically.  The  type  of  factors  referred  to  include  FORTRAN  do-loop 
overhead,  computer  data  access  times,  etc.  To  this  end,  two  digital 
computer  programs  were  written  in  FORTRAN  IV  for  the  CDC-6400  computer 
which  were  completely  Identical  except  for  the  solution  algorithm.  In 
one  program  the  array  algebra  algorithm  was  implemented,  and  in  the  other 
the  standard  least-squares  normal  equations  were  utilized. 

The  first  program,  array  algebra  solution  (AASOL) , performs  the 
following  operations: 

1.  Generates  an  elevation  grid  (Z,)  which  has  an  exact  least- 

A 

squares  approximation. 

2.  Generates  the  matrix  (X)  of  appendix  C and  finds  the  expression 

T -1  T 
(X  X)  ^X 

3.  Solves  the  equation 

= (x'’^x)"^x'^z^y(yTy)“1 

Note  that  since,  in  this  application  X = Y,  then 

T 

Y(Y^Y)~^  = t(X^X)"^x'^ 

and 

T -1 

Y(Y  Y)  need  not  be  calculated. 
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4.  Evaluates  the  matrix  of  coefficients  and  produces  a 

set  of  derived  elevation  data.  This  data  is  then  compared  to 
and  a sum-squared  error  computed. 

5.  Measures  the  CP  times  of  the  one-time  run  of  step  (2)  and  the 
repetitive  run  of  step  (3),  and  displays  the  total  computation  time 
needed  to  model  a complete  1:50,000  map  sheet. 

Successful  runs  were  achieved  for  model  equations  of  4,  9,  16,  25, 

36,  49,  64,  81,  and  100  terms,  with  computation  times  exceeding  the 
multiplication  and  addition  times  (as  predicted  in  the  computational 
analysis)  by  a factor  of  approximately  2.7.  The  actual  CP  execution 
times  ran  from  166  seconds  for  the  4-term  model  polynomial  to  806 
seconds  for  the  100-term  model. 

The  conventional  Least-Squares  Solution  test  routine,  LSSOL,  performed 
functions  equivalent  to  the  Array  Algebra  Solution  test  routine. 

Successful  runs,  however,  could  be  achieved  only  for  model  polynomials  of 
4,  9,  and  16  terms.  Any  attempt  to  dimension  the  routine  to  accept  a 
larger  model  equation  exceeded  the  storage  capablty  of  the  CDC-6400. 
Computation  times  exceeded  the  predicted  multiplication  and  addition 
times  by  a factor  of  2.3.  Actual  CP  execution  time  ran  from  218  seconds 
for  the  4-term  model  to  941  seconds  for  the  16-term  model.  The  comparison 
of  the  computation  times  in  AASOL  and  LSSOL  are  presented  in  figure  5. 

The  savings  ratio  actually  achieved  for  the  most  efficient  4- 
coefflclent  model  case  of  our  terrain  modeling  procedure  was  1.3  to  1, 
more  than  50  percent  lower  than  the  theoretically  projected  ratio  of 
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Minutes  of  CDC  6400 
computer  time  per  map 
sheet  required  to 
calculate  initial 
polynomials. 


1 1 1 1 

4 16  36  6- 

Number  of  terms  in  model  polynomial 


Figure  5:  Empirical  comparison  of  efficiency; 
Least-Sauares  vs.  Arrav  Aleebra 
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1.75  to  1.  This  difference  is  caused  by  the  empirically  determined  higher 
implementation  overhead  associated  with  the  array  algebra  algorithm — 

2.3  times  the  least-squares  projection  vs.  2.7  times  the  array  algebra 
projection.  The  higher  implementation  overhead  of  array  algebra  appears 
to  be  due  to  the  much  larger  number  of  do-loops  required.  The  major 
result,  however,  is  that  the  array  algebra  approximation  algorithm  is 
more  efficient  in  all  cases. 

Examining  the  possible  savings  at  a production  center  if  the  array 
algebra  algorithm  were  used  in  the  most  efficient  case; 

1.  52  seconds  were  saved  per  map  sheet 

2.  500  map  sheets  are  produced  per  year^^ 

3.  $160  per  hour  is  the  current  cost  of  computer  processor  time 
at  a production  center.  Implementation  of  the  array  algebra  algorithm 
will  therefore  result  in  an  annual  savings  of  $1,155.56.  Since  our  best 
estimates  indicate  that  array  algebra  will  require  at  least  three  man-months 
to  Implement  in  CONSAC  II,  such  an  implementation  cannot  be  considered  cost- 
effective  when  viewed  only  in  terms  of  its  computational  efficiency.  There 
is  an  effort  presently  under  way,  however,  to  modify  CONSAC  II  to  utilize  a 
9-term  polynomial  model  in  order  to  investigate  possible  improvement  of 
contour  fitting  and  expression.  Array  algebra  can  be  economically  included 

in  this  effort  and,  if  used  in  production,  will  result  in  considerably  greater 
savings  than  if  used  with  a 4-term  model. 
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Conclusions. 


An  investigation  of  the  applicability  of  Rauhala’s  array  algebra  to 
the  ETL  terrain  modeling  algorithm  has  been  completed.  Our  theoretical 
analysis  of  the  two-dimensional  bilinear  polynomial  case  has  verified  its 
equivalence  with  the  least-squares  normal  equations  for  orthonormal 
observation  grids  and  provides  for  a more  traditional  and  illuminating 
derivation.  The  comparison  of  the  computational  requirements  paper  study 
and  the  empirical  analysis  has  shown  that: 

1.  the  array  algebra  algorithm  has  significantly  more  implementational 
overhead  than  the  least-squares  algorithm. 

2.  the  array  algebra  algorithm  is  more  efficient  than  the  least- 
squares  algorithm  in  all  cases  pertinent  to  our  requirements. 

3.  for  our  currently  employed,  most  efficient  terrain  modeling  case, 
array  algebra  has  an  improvement  ratio  of  1.3  to  1,  amounting  to  a 
savings  of  only  about  $1,000  per  fiscal  year  in  a production  environment. 
Based  on  our  analyses,  implementation  of  the  array  algebra  algorithm  in 

the  specific  case  investigated  is  not  as  significant  as  had  been  anticipated. 

Our  current  plans  are  to  proceed  with  the  array  algebra  implementation 
at  the  same  time  that  we  are  expanding  the  terrain  modeling  software's 
capability  under  a related  R&D  effort Modification  of  our  terrain 
modeling  software  to  the  array  algebra  algorithm  can  be  economically 
included  in  planned  modifications  under  that  work  effort  giving  us  the 
opportunity  to  gain  further  experience  with  this  powerful  algorithm, 
which  we  feel  can  be  profitably  employed  whenever  approximation  of 
observations  on  an  orthonormal  grid  is  required. 
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APPENDIX  A.  THEORETICAL  ANALYSIS  OF 
ARRAY  ALGEBRA  AS  APPLIED  TO  ETL’S 
TERRAIN  MODELING  ALGORITHMS 


Preliminaries . The  problem  statement  proceeds  as  follows:  It  is  desired  to 
approximate  the  q by  p set  of  observations, 


at  , Y.  for  i = l,  2,  3,  ...  , q 

j=l,  2,  3,  ...  , P 


where  the  X^  and  Y^  are  the  independent  variable  locations  of  the  measured 
observations,  (It  is  important  to  note  that  the  lie  on  an  ortho-normal  grid 

in  the  plane  of  the  independent  variables  X,Y  by  construction.')  with  the  model 
equation. 


W 


ij 


T 

Z 

r*0 


U 

Z C x5  Yf 
rs  i J 


Where  the  C^g  are  the  unknown  approximating  function’s  coefficients,  where  (T+1) 
times  (j+1)  is  less  than  Q times  P.  It  is  now  evident  that  the  model  equation  can 
be  written  in  matrix  form  (using  the  bilinear  form)  as: 


[1,  X^, 


1x1  lx  (T+1) 


So’  ^01’  •••  ’ ^ou 

1 

^10’  ^11’  •••  ’ S.U 

,.2 

• • • 

• • • 

• 

Cto  * « • • • » 

• • 

1 

(T+1)  X (U+1)  (I'+l)  X 1 
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where  the  size  of  the  matrices  is  shown  underneath  each  matrix. 


Further,  it  is  possible  to  define  observational  matrices  which  include  all  of 
the  observations,  e.g. , 


s 

w « 

'Wii  , 

^12  * • • * ’ 

W21  , 

^22  * • “ • * 

'^2p 

• 9 

• 9 9 

• 

W , 
<11 

Wq2  » . • . > 

and  then  the  entire  set  of  model,  observation  equations  can  be  written  in  matrix 
form  as: 


W = X C y"^ 

qxp  qx(T+l)  (W-1)  x p 
(T+1)  x (Uhl) 

where 
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and 


1 . . Yt  , 


1 , ^2  , , 


1 , Y , Y^ 
P P 


Y^ 


One  other  preliminary  definition  necessary  in  the  following  sections  is  the 
differentiation  of  an  arbitrary  scalar  function,  with  respect  to  an  arbitrary 
rectangular  matrix,  D , where 

®11  » ®12  » *^13  » • • • » 


D - 


°21  ’ ^22  . ^23  , . . . , D23 


°ml  ♦ °m2  ♦ °m3 °ml 

as  is  customary  the  matrix  of  partlals  is  denoted  as: 


li 

3D 


li 

li 

36 

3Dii  , 

^^12  • • ’ * • 

3Dil 

li 

li 

li 

3D22^  , 

3D22  » • • • » 

3D21 

33 

• • t 

li 

36 

3Dml  • 

^°m2  • • • • » 
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4.  Derivation  of  Array  Algebra  Solution;  There  are  an  infinite  number  of  solutions 
to  the  approximation  problem  as  stated  in  section  three.  A unique  solution  can 


only  be  achieved  if  additional  conditions  are  specified.  The  condition  which  has 
been  determined  (by  the  study  reported  upon  herein)  to  result  in  Rauhala’s  array 
algebra  solution  is  as  follows: 

Determine  the  coefficients  of  the  model  equation,  C^g  , which  minimize  the  scalar 
function,  (fi  , defined  as: 

4.  - II  w - X c y11 


where  the  symbol  ll'-ll  is  the  squared  matrix  norm,  which  is  defined  for  an 
arbitrary  rectangular  mxl  matrix,  D , as 


m 

Z 

1=1 


The  careful  reader  will  note  that  this  is  a simple  reformulation  of  the  scalar 
objective  function  utilized  to  derive  the  standard  least-squares  solution) 

This  reformulation  of  the  least-squares  objective  function  simply  utilized  the 
bilinear  form  characteristic  of  and  required  by  array  algebra  and  the  well-known 
concept  of  matrix  norms. 
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Proceeding  with  the  derivation,  the  objective  function  is  minimized  bv  differen- 
tiating it  with  respect  to  the  unknown  model  coefficients,  C^g  and  setting  the 
result  equal  to  zero: 

T U 

M.  = i ^ ^ 

k=0  1=0 


where : 


then: 


Wki  - [1.  \ . 


C 1 


T u ae^i  _ T u ^ s 

H E 201,1  ac.,  ■ E r 29  , (-Xk  Yi) 

i,,o  1=0 


k-0  1=0 


T U 

-2  E E Xk  01,1  Yj 

k=0  1=0 


-2  (Xi  , X2 xb  [W-X  Cy'^] 


1^  = -2  x'^(W-X  CY^)Y 
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and 


-2  X^[W-X  C Y^]Y  = 0 

multiplying  the  terms  out,  dividing  by  -2  and  rearranging  gives: 

x'^X  C y’^Y  = X^  W Y 

T — 1 T 

postulating  the  existence  of  (X  X)  and  (Y  Y)  , and  multiplying  them  on 
the  left  and  right  respectively  gives: 

C = (X'^X)"^  x'^  W Y (Y^Y)”^ 


which  is  exactly  Rauhala's  array  algebra  solutioni  It  is  now  obvious  that  the 
array  algebra  solution  must  give  exactly  the  same  numerical  results  as  the  least- 
squares  solution  for  those  approximation  problems  which  have  their  observations 
on  an  orthonormal  gird — they  must,  they  are  both  the  uniquci  solution  of  over- 
determined systems  subject  to  exactly  the  same  conditions! 

5.  Approximation  with  Orthonormal  Observational  Weighting;  The  capability  for 
nearly  arbitrary  weighting  of  the  observations  is  a requirement  for  many  of  the 
MC&G  approximation  applications.  After  considerable  analysis,  it  was  determined 
that  completely  arbitrary  weighting  of  the  observations  was  impossible  if  the 


bilinear  form  and  array  algebra  solution  were  used.  A compromise,  orthonormal 
weighting,  however,  was  successfully  found  which  meets  the  observational  weighting 
requirements  of  our  applications.  By  orthonormal  Weighting,  It  is  meant  that  a 
single  unique  weight  can  be  associated  with  each  x and  each  y coordinate 
specifying  the  observation's  ortho-normal  grid.  The  weight  corresponding  to  any 
single  observation  Is  the  product  of  Its  x-welght  and  y-welght.  This  scheme  Is 
Implemented  as  follows:  The  objective  function  is  modified  to: 

^ = II  “x  ^ “yll 


where 
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and  the  weighting  associated  with  the  i.j-th  observation  is  given  by 

OJ  = CO  CO 

ij  y. 

This  scheme,  while  not  completely  arbitrary,  does  allow  for  great  flexibility  in 
weighting,  with  the  observation  weighting  over  the  rectangular  area  being  the 
product  of  the  x and  y variances.  The  derivation  for  the  weighted  solution  now 
follows  as: 

objective  function:  4>  » ||njj  (W  - X C Y^)  fiyjl 

minimize  by  setting  partials  with  respect  to  unknowns  to  zero; 

=•  -2  x’’  n [W-X  C y'^]  n Y = 0 

3c  X . y 


rearranging: 


x’^  Wfiy  Y - x'^fix  X C y"^  fiy  Y =*  0 


postulating  existence  of  inverses: 

c = (x”^  ^x  x)"^  x'’^  w Y(Y'^nyY) 


This  solution,  as  in  the  unweighted  case  of  section  four,  must  produce  exactly 
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the  same  numerical  results  as  the  standard  least-squares  solution  with  the  same 
weighting — because  both  techniques  produce  the  unique  solution  based  on  exactly 
the  same  mathematical  critera. 

6.  Constrained  Array  Algebra;  As  mentioned  earlier,  it  is  often  necessary  to 
rigorously  constrain  the  approximating  function's  value  at  more  precisely  known 
points.  Utilizing  the  matrix  norm/bilinear  form  objective  function  and  the  cal- 
culus of  variations  with  Lagrange  multipliers,  observation  constraints  can  be 
rigorously  incorporated  into  the  array  algebra  solution. 

The  derivation  of  the  constrained,  weighted  array  algebra  solution  proceeds 
follows : 

minimize  the  objective  function;  <)>  , 

<t>  (W  - X C y’’^) fly- 

subject  to  the  N constraints: 

'f'k  * *ck  ^ ^ck  k = 1,  2,  3 N 
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where: 


X 


ck 


[1. 


ck’ 


Xckl 


Where  is  the  X location  of  the 

k-th  observation  constraint,  . 


- [1»  ^ck’  ^ck»  • • • » ^ck^ 

Where  is  the  Y location  of  the 

k-th  observation  constraint,  . 


NOTE:  The  location  of  the  constraints  are  not  restricted  to  the  orthonormal 
observation  grid  values! 

utilizing  the  modified  objective  function,  (|) ' : 

N 

=||“x  (W  - X C YT)Qy||+  Z - X^j^  C Y^j^) 

k”  1 

where  are  the  N Lagrange  multipliers. 

Minimizing  this  function  gives: 

t N 

g - -2  X^n^[W  - X C Y^lfiy  Y - X^^  Y^^ 

solving  for  C ; 

T T N „ 

-2  x^n^  w n Y + 2 X X c Y ^ ^ \ » o 

^ k=l 
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rearranging: 


X”’'^  fix  * C Y 


N 

fi^  W fi„  Y + E \ 
^ k-1  2 


assuming  the  Inverses  exist: 


-1 


N 


c - (x^  X)  [X  (y"  fiy  y) 


-1 


and  to  evaluate  the  Lagrange  multipliers: 

♦k  ■ *ck  k - 1.  2.  3.  ...  H 


substituting: 


♦k 


X (x"^  fi^  X)"^  [x’^  fi„  W fi^  Y + E ^ x’’^  Y 1 (Y^  fi^  Y)~^  Y^ 
ck  * * y cj  cj  y ck 


Since,  by  design,  there  are  N constraint  equations  and  N*-  Lagrange  multipliers, 
this  set  of  equations  represents  a conslstant,  determined  set  that  can  be  solved 
for  the  unknowns. 

The  closed-form  matrix  equation  for  the  Lagrange  multipliers  In  the  more  general 
case  N>1  has  not  been  found  as  yet. 
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Consider  the  case  when  only  one  constraint  is  present;  then,  the  constraint 


equation  becomes 


“1  T 

iP  = (x'^  X)"^  [x"^  w fiy  Y]  (y'’^  Y)  Y^ 

+ A X^  (X'^  xl  Y^  (Y^  fiy  Y)’^  Y^ 


or  rearranging: 


ij;  - X (x'^  n x)"^  [x’^  n w n y]  (y"^  n y)"^  y 
c X X y y c 

x^  (xT  X)-1  x'J  Y^  (yT  n y)"^  y^ 
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APPENDIX  LEAST-SQUARES  POLYNOMIAL  FITTING 


THE  DERIVATION  OF  = (a'^A)  “^A^Z^ 


n 

Given  a set  of  n elevation  observations,  I , suppose  that  one 
wishes  to  fit  the  polynomial  i=l 


fCx^.y^)  = Cq  + + C2y.  + C3X.y^ 


as  closely  as  possible  to  these  observations.  An  accepted  "best"  fit 
occurs  when  the  sum  of  the  squares  of  the  vertical  distances  from  the 
observations  to  the  described  surface  is  minimized.  The  distance  from  an 
elevation  point  to  the  described  surface  may  be  defined  as 


d = - f(x^,y^) 

= (z^  - f(x^,y^))2 

The  set  of  n distance-squares  is  therefore 

n 2 " 2 

S = Z di  = Z (='•1  ■ Cq  - c^x.  - cyy,  - C3X.y.) 
i=l  i=l 

f 

S will  be  minimized  when  its  derivative  is  equal  to  zero. 
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* <SS  = -2E(z^  ~ ^0  ~ “ ^2^1  “ ° ° 


6c 


6S 


6c, 


6S 


6c_ 


6S 


6c, 


-2E[(z^  ~ <^0  " " ^^2^1  “ C35'iyi)  (’'i)  1 = 0 


-21  [(z^  ■ " *^1^1  ■ ‘■■2^1  " c^x^y^Xy^)]  =0 


= -2Et(z^  “ *^0  " ■ *^2^1  ~ C3X^y^)(x^y^)]  = 0 


> 


(Bl) 


CQn  + c^  E X.  + C2  E y.  + C3  E x^y^  = Ez^ 

CqExj^  -f-  Cj^  Ex^  + c^  ■*■  *^3  ^^fyi  ~ ^x^z^ 

Cq  E yi  + c^  E x^  y^  -H  C2  E y^  + c^  E x^  y^  = E y^z^ 

Cq  Z x^y^  + c^  E x2y^  -H  C2  Z x y^  -H  C3  E xjy^  = e x^y  .z^ 


> 


(B2) 


*Hereafter  the  symbol  E indicates 


n 

E 

1=1 
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Define  A as  an  n by  n matrix  of  polynomial  terms,  as  nn  m by  1 

matrix  of  coefficients,  and  as  an  n by  1 matrix  of  observations. 

i 

Then, 
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The  left  side  of  equation  (b2)  noy  be  written  as 


which  in  matrix  notation  is 


aX 

The  right  side  of  equation  (B2)  can  now  bo  written  as 
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n * « _ 


I 


z(0,0)  = 100 
2(0,1)  = 110 
z(0,2)  = 112 
2(1,0)  = 118 
z(l,l)  = 108 
z(l,2)  » 120 
z(2,0)  = 106 
z(2,l)  = 104 
z(2,2)  = 114 
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The  model  equation  is  chosen  as 


f(x,y)  = Cq  + CjX  + c^y  + c^xy 


9 

9 

9 

9 

= 9 

15 

9 

15 

9 

9 

15 

15 

9 

15 

15 
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10  0 0 
10  10 
10  2 0 
110  0 
A = 1 1 1 1 

112  2 
12  0 0 
12  12 
1 2 2 A 
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9 
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25 

10 

-5 

10 

A 

-2 

-5 
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-15 

-6 

3 
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0 

0 

15 

6 

-15 

0 

15 

-6 

0 

6 

3 

0 

9 

0 

-9 

0 

0 

0 

-9 

0 
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Cq  ■ 105.222  ^ C2  ■■  4 2/3  • -1 

s - 105.222  * 1.333X  4.667y  - xy 

!foclc«  that  this  rasult  is  Idsntlesl  with  that  of  the  array  algebra 
luwerical  exaaple  in  appendix  0« 
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APPENDIX  C.  ARRAY -ALUEllRA  POLYNOMIAL  FITTING 


THE  DERIVATION  OF  C = (x'^X)"^x'''z*Y(Y^Y)"^ 

A A 

When  the  set  of  elevation  observations 


m-1 


is  ordered  in  an  orthonoriral  grid,  then  It  may  be  written  as 


(a-l),(b-l) 

£ Z4 


l,J-0 


ij 


with  ab=m  terms  in  the  model  equation.  If  the  model  equation  is  chosen 
as  the  four-term  polynomial 


- (Xi-Yi)  = CqO  ‘^10*i  ‘^OlY^  + ciix^y^ 

then. 
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In  the  grid  are 


and,  given  n*ef  observations,  then  all 


= Z*  - XC 
A A 


where 


It  is,  of  course,  possible  to  generalize  this  derivation  by  selecting  a 
general  model  polynomial. 
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Thus,  choosing 


. .+c 


(a-l)(b-l)-  i 


.(a-l)„(b-l) 


j 


then. 


(a-1) .(b-1) 


^±i- 


2 3 

"l  *i- 


— 

o 

o 

o 

^01 

^^02 

' • ‘^O(b-l) 

1 

‘^lO 

"ll 

Ci2 

* • ‘^Kb-D 

yj 

2 

o 

CM 

O 

^^21 

^^22 

• • ‘^2(b-l) 

• 

• 

• 

• 

• 

« 

• 

• 

• 

• 

• 

• t 

• 

„(b-n 

‘'(a-l)O 

‘'(a-l)l 

''(a-l)2  • 

' ■ ‘^(a-l)(b-i: 

''j 

46 


and  all  In  the  grid  are 


1 .. 


1 XjX2 


1 Xj  Xj 


x(*~l)l 


.(•-1) 


x(«-l) 


1 X X*  ...  x^-lH 
A e 


00 


10 


■20 


■'Ol 


21 


'02 


"12 


'22 


®0(b-l) 

*l(b-l) 

®2(b-l) 


'(•-1)0  *'(,-1)1  '(•-1)2”*  '(a-l)(b-X^ 


y\ 


2 

>3 


...  1 


2 

yf 


y(b-l)  y(b-l)  y(b-l)... 


Za  - XCaY 


^ere  X la  an  e by  a matrix  of  x-dlrection  parameters,  Y is  a f by  b 
natrix  of  y-directlon  parameters,  Ca  la  an  a by  b matrix  of 
coefficients,  and  Za  is  an  e by  f matrix  of  elevation  observations. 

Solving  for  C. 

A 

h ■ 

Z.Y  « X C y’’y 
A A 

x'^Z.Y  - (x’^X)C  (Y^Y) 

“ A 

(X^X)"^  x'^’ZaY  - (x'^X)"^(X^X)Ca(y’’y)  - Ca(y'’^Y) 
(X^X)"^x’^ZaY(y'^Y)''^  - Ca(Y^Y)(y'^Y)"^  - 
Ca  - (x^x)"^x^ZaY(y’'y)"^ 


Note  that  when  X * Y,  a special  computational  case  exists: 
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y(y’^y)‘^  = X(X^X)"^ 


= [^lx('^x)"^]‘^j^ 
= [^[(x'^x)‘^]Vj 


T —1 

and,  since  (X  X)  is  symmetric. 


(x'^X)"^]Vj 


T -1  T, 

= [(x'x)  ^X  ] 


therefore,  when  X = Y , 

= [(X^X)"^  x"^]  [(X^X)"^  x'^]^ 

T 

and  Y(Y  Y)  need  not  be  calculated. 

A NUMERICAL  EXAMPLE  OF  THE  ARRAY  ALGEBRA  TECHNIQUE.  The  same  observational 
values  used  in  the  numerical  example  in  appendix  A will  also  apply  to  this 
example. 

The  model  equation  z = c^^  + Cqj^x  + + Cj^^xy  can  be  written 


z - (l,x) 


00  *^01 


10  *^11 


therefore. 


‘‘OO  *"01 


^10  ‘^ll 


(x^x)"^  x'^z.y(y'’^y)"^ 
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T 1 T -1  T ^ 

since  Y = X,  Y(Y  Y)“^  = [ (X  X)  X ] 


_1_  5-3 

Y(Y  Y)"^  =6  2 0 

-1  3 

1 5 2 -1  100  110  112  '_1_  5 -3 

C = 6 -3  0 3 118  108  120  6 2 0 

^ 106  104  114  -1  3 

• j J L « 

630  662 

= 36  18  -18 

^ 3788  168 

■ 36  48  -36 
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105.222  4 2/3 

= 4/3  -1 

Z = 105.222  + 1.333X  + 4.667y  - xy 


Notice  that  this  result  is  identical  with  that  of  the  least-squares 
numerical  example  in  appendix  B. 
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APPENDIX  D . THE  NUMBER  OF  MULTIPLICATIONS  AND  ADDITIONS 
REQUIRED  TO  FIT  A POLYNOMIAL  OF  m TERMS 
TO  DATA  SETS  OF  n OBSERVATIONS 


REQUIRED  MULTIPLICATIONS  AND  ADDITIONS  USING  THE  METHOD  OF  LEAST-SQUARES. 
Given  that 

- (a’‘a)~^a’'  Zl 


tfhere 


A is  an  nxm  matrix,  Z la  an  nxl  matrix. 

la 


Then 


T 

L ■ A A requires 


nZl  multiplications  and 

1-1 


whic^  simplify  to 


— nm(i»4-l)  multiplications  and 


m 


(n-1)  El  additions 
1-1 


m(n-l)'(mfl)  additions. 


P - (L)”^  requires 


1 2 

— m(4m  +3m-l)  multiplications  (see  reference  13) 
6 


and 


1 2 

•g-  m(4m  -3m-l)  additions 


(see  reference  13) 
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T 

Q = PA  requires 

m n multiplications  and  ran(m-l)  additions. 

R = (a'^A)~^a'^  requires 

J.  m(4m^  + 3m  + 9mn  + 3n  - 1)  multiplications 

6 

and 


1 9 

_ m(4m  - 6m  + 9mn  - 3n  - 4)  additions. 

6 


T -1  T 

In  the  polynomial  fitting  process,  it  is  necessary  to  compute  (A  A)  A 
only  once.  The  multiplication  of  (A'^A)"^A'^  to  Z^,  however,  must  be 
repeated  as  many  times  as  necessary  to  compute  all  of  the  polynomials 
needed  to  describe  the  terrain. 

This  process  requires 

(mn)  (//  of  reps.)  multiplications*,  and 

m(n-l)(#  of  reps.)  additions. 

REQUIRED  MULTIPLICATIONS  AND  ADDITIONS  USING  ARRAY  ALGEBRA. 

= (X'^X)"^X'^  Y(Y^Y)"^ 

where 


X 

is  an 

e 

X a 

matrix 

Y 

is  an 

f 

X b 

matrix 

c* 

ah  = m (terms  in  the  model  equation) 


5Z 


T 

X X requires 


4.  ac(a  + 

7 


1) 


and 


- a(e  - l)(a  + 1) 
2 


multiplications 


additions. 


(X^X)”^  requires 


i.  a(4a^  + 3a  - 1) 
6 


and 


multiplications  (see  reference  13) 


3a 


1) 


additions  (see  reference  13) 


(X^)"^  X^  requires 

a^e  multiplications 

and  ‘ 

ae(a  - 1)  additions. 

(X^X)”^X^  therefore  requires 

1 a(4a^  + 3a  + 9ae  + 3e  1)  multiplications 
6 


and 
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additions. 


i-  a(Aa^  - 6a  + 9ae  ~ 3e  - 4) 
6 


Y(Y^Y)“^  likewise  requires 


— b(Ab^  + 3b  + 9bf  + 3f  - 1)  multiplications 
6 


- b(4b^  - 6b  + 9bf  - 3f  - 4)  additions. 


T -It 

The  repetitive  computations  involve  the  multiplication  of  (X  X)  X 
to  Z and  the  multiplication  of  (X^X)~^X^Z  to  Y(Y^Y)~1.  These  steps 
require 


p(an  + mf) 


multiplications 


- p(an  - af  + raf  - m) 


additions 


where  p is  equal  to  the  number  of  polynomial  fits  per  map  sheet. 
Where  X = Y, 


T -1  T -1 

Y(Y  Y)  ^ = X(X  X)  ^ 


f[X(x'^x)-‘ 


which,  since  (X^X)~^  is  symmetric,  is  equal  to  [(X^X)"'x'*^] 


T 1 

therefore,  Y(Y  Y)“^  need  not  be  ca] ciil.Tted . 
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APPENDIX  E.  POLYNOMIAL  COVERAGE  AND  THE 
NUMBER  OF  POLYNOMIALS  PER  MAPSHEET 


Reference  two  presents  empirical  evidence  that  the  best  approximation  of  automatic 
compilation  elevation  data  for  a 1:50,000  output  scale  and  a density  of  approxi- 
mately one  million  points  per  mapsheet  Is  achieved  for  a roving,  locally  valid, 

2 

four  coefficient  polynomial  model  equation  when  It  Is  calculated  from  169  = 13 
data  points  using  46Z  overlap.  For  the  purpose  of  studying  larger  model  equattoas 
It  was  assumed  that  each  coefficient  added  would  allow  for  a proportional  Increase 
In  area  covered,  e.g.,  the  number  of  terms  In  the  model  equation,  m , would  be 
related  to  the  nundier  of  observations  used  to  compute  It  by  the  formula: 

n - (INTEGER  V»2.25m)^ 

which  Is  derived  directly  from  the  case  In  point. 


Using  the  46Z  overlap  and  a typical  916  x 1112  “ 10^  data  grid  for  a 1:50,000 
mapsheet,  the  number  of  polynomials,  p , needed  per  map  sheet  Is  given  by: 


P 


[INTEGER  ROUNDED  UP 


-ilL.  ][ INTEGER  ROUNbED  UP 
>s(e+l)" 


1112  , 
%(f+l) ^ 


where  e ” f * VS-  for  unbiased  modeling.  The  denominators  compute  the  point 
spaclngs  between  neighboring  local  polynomial's  centers.  The  two  factors  used 
then  compute  the  number  of  rows  and  columns  of  polynomials,  respectively.  Note 
that  by  successive  use  of  these  formulae,  the  effect  of  larger  polynomial  model 
equations  resulting  In  fewer  polynomials  per  sheet  Is  automatically  accounted  for. 
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